AD-A090  170  SKF  TECHNOLOGY  SERVICES  KING  OF  PRUSSIA  PA  F/0  V 

TWO-WAY  ANALYSIS  OF  VARIANCE  FOR  WE I BULL  POPULATIONS. (U> 

AUG  80  J  I  MCCOOL  F49620-79-C-0035 

UNCLASSIFIED  SKF-AT80D048  AFOSR-TR-80-0999  NL 


wiur.LtoWU  AD  A  0  9  0 1 7 


FINAL  SCIENTIFIC  REPORT 


TWO-WAY  ANALYSIS  OF  VARIANCE  FOR  WEIBULL  POPULATIONS 


AUGUST  1980 


J.  I.  MC  COOL 


SKF  REPORT  NO.:  AT80D048 


p  ^  \ 


SUBMITTED  TO: 

DIRECTOR,  MATH.  &  INFO.  SCIENCES 
AIR  FORCE  OFFICE  OF  SCIENTIFIC  RESEARCH 
ATTN:  NM,  BLDG.  410 
BOLLING  AFB ,  D.C.  20332 

CONTRACT  NO.  F49620-5-C-0035 


SUBMITTED  BY: 

SKF  TECHNOLOGY  SERVICES 
~  SKF  INDUSTRIES,  INC. 
KING  OF  PRUSSIA,  PA. 


Approved  fop  public  release  | 
distribution  ualiaited. 


80  10  6  083 


secoaiTY  class*  pi  cation  op  this  page  (Wh«« o atm  Entatad) 


I  ^REPORT  DOCUMENTATION  PAGE 


SK 


m 

ifa 


OSRyTR-  8 


TBESBSk 


2.  GOVT  ACCESSION  NO.I 


Ac#  oil' 


f.  TITLE  (and  Subtltla) 


5.  TYPE  OF  RERO' 


Two-Way  Analysis  of  Variance  for  I 
Weibull  Populations  »  j — — ' 


Scientific,  Final 
12/30/78  -  7/31/80 


•.  PERFORMING  ORG.  REPORT  NUMBER 


a.  contract  on  grant  number/ •) 

F  4  9  6  2/-  7  9  -  C  -jlj/l 


9.  PERFORMING  ORGANIZATION  NAME  AND  AOORESS 

SKF  Industries,  Inc. 

1100  Pirst  Avenue  / 

King  of  Prussia,  Pa.  19406  ' 


11.  controlling  OFFICE  NAME  ANO  AOORESS  A 

Air  Force  Office  of  Scientific  Reseat 
Bolling  AFB ,  D.  C.  20332 


|  ;  Au 

.  _  i  1 ' 


Ml  1  23 


IS.  SECURITY  CLASS,  (el  thlt  report; 


Unclassified 


15*.  OECLASSI  Fl CATION/ DOWNGRADING 
SCHEOULE 


IS.  DISTRIBUTION  STATEMENT  (of  (hit  Report) 


Approved  for  Public  Release:  Distribution  Unlimited 


17.  DISTRIBUTION  STATEMENT  (o (  tba  abatract  aatarad  In  Block  20,  it  dllfarant  from  Raport) 


19.  KEY  WORDS  fContlnua  on  ravaraa  aid*  If  nacaaaary  and  idantlfy  by  block  numbar) 


Weibull  distribution,  analysis  of  variance,  maximum  likelihood 
estimation ,  single  classification,  regression  analysis,  location 
parameter,  grouped  data 


20  A  B&T R Afc T  (Contlnua  on  ravaraa  aids  It  nmcaaaary  and  Idantlfy  by  block  numbar) 


Expressions  are  given  for  the  maximum  likelihood  estimators 
of  the  Weibull  shape  parameter  using  data  from  a  two-way  class¬ 
ification  experiment  under  five  hypothesized  relationships  for 
the  scale  parameter  in  each  cell,  viz,  the  scale  parameter  (1)  is 
constant  (2)  has  a  multiplicative  row  effect  (3)  has  a  multipli¬ 
cative  column  effect  (4)  has  a  multiplicative  row  and  column  j 
effect,  and  (5)  has  multiplicative  row  and  column  and  inter- \  I 


Unclassified 


3BL 


9 


LASSIFICATION  OF  THIS  PA<je<T*»i*«  Enttfd) 


20. 


action  effects.  Procedures  for  analyzing  a  two-way  lay¬ 
out  without  interaction  (a  randomized  block  design)  are 
discussed  and  illustrated. 


Techniques  are  given  for  estimation  of,  and  infer¬ 
ence  on,  the  location  parameter  of  the  three-parameter 
Weibull  distribution  and  a  military  application  of  the 
technique  is  given 

Software  improvements  are  cited  in  the  Monte  Carlo 
simulation  computer  program  for  the  computation  of  the 
distributions  of  pivotal  functions  required  for  analyzing 
Weibull  regression  experiments. 

In  these  experiments,  the^response  variable  follows 
a  two-parameter  Weibull  distribution  with  a  scale  para¬ 
meter  that  varies  inversely  with  a  power  of  an  external 
variable  termed  a  "stress."  With  the  software  improve¬ 
ments,  pivotal  functions  are  determined  which  permit  (1) 
testing  the  adequacy  of  the  power  function  fit,  (2)  set¬ 
ting  confidence  limits  on  the  shape  parameter  and  stress- 
life  exponent  and  (3)  setting  confidence  limits  on  a  spe¬ 
cified  percentile  of  the  distribution  at  seven  prescribed 
"stress"  levels.  Tables  have  been  generated  for  the 
analysis  of  experiments  conducted  at  two  and  three  stress 
levels,  with  sample  sizes  ranging  from  5  to  20. 

Methodology  for  the  analysis  of  a  single  classifi¬ 
cation  experiment  with  Weibull  distributed  response  is 
described  and  the  contents  of  an  applications  oriented 
paper  is  described.  Other  applications  of  the  methodol¬ 
ogy  are  cited. 

Software  is  described  for  simulation  of  the  distri¬ 
bution  of  the  maximum  likelihood  estimates  of  the  Weibull 
parameters  that  result  when  the  data  are  grouped.  Results 
of  varying  interval  width  and  sample  size  are  given. 


(Uk£JJim£l£R 


f.  >i.> safe.*;.*  > 


PROPRIETARY  NOTICE 

11 

Furnished  under  U.S.  Government  Contract  No.  F4962O-|0-C-OC35 . 
Shall  not  be  either  released  outside  the  Government,  or  used 
duplicated,  or  disclosed  in  whole  or  in  part  for  manufacture  or 
procurement,  without  the  written  permission  of  the  Director, 
Math.  &  Info.  Sciences,  Air  Force  Office  of  Scientific  Research, 
Bolling  AFB,  D.C. 


(X(  ** 


<3? 


-V,  X 


AIR  FORCE  OFFICE  OF  SCIENTIFIC  wcno/iu  _ 
NOTICE  OF  TRANSMITTAL  TO  DDC  £ARCH  (ArSCi 
This  technical  report  ►vv*  k ...  < 
approved  t  or  put>l  i o  rc  lea'-s  rim  j™,8*94  and  13 
Distribution  is  unlimited'*  190-18  (7b), 

A.  D.  BLOSE 

Teohnloal  Information  Offioer 


SKr  TtCHNCLCX^v  SERVICES 


"A?  RUSSIA  PA 


onui  Km* 


FINAL  SCIENTIFIC  REPORT 


TWO-WAY  ANALYSIS  OF  VARIANCE  FOR  WEIBULL  POPULATIONS 


AUGUST  1980 


J.  I.  MC  COOL 


PREPARED 


APPROVED: 


RELEASED: 


SKF  REPORT  NO. :  AT80D048 


SUBMITTED  TO: 

DIRECTOR,  MATH.  &  INFO.  SCIENCES 
AIR  FORCE  OFFICE  OF  SCIENTIFIC  RESEARCH 
ATTN:  NM,  BLDG.  410 
BOLLING  AFB ,  D.C.  20332 
CONTRACT  NO.  F49620-*P-C-0C35 

n 


SUBMITTED  BY: 


)KF  TECHNOLOGY  SEPVCES 


■;f  :  ’  Z'  c'c‘ 


AT80D048 


TABLE  OF  CONTENTS 


1.0  SINGLE  CLASSIFICATION  EXPERIMENTS  1 

2.0  LOCATION  PARAMETER  ESTIMATION  AND  INFERENCE  2 

3.0  WEIBULL  REGRESSION  OR  ACCELERATED  TEST 

ANALYSIS  5 

4.0  TWO-WAY  CLASSIFICATION  EXPERIMENTS  8 

5.0  GROUPED  DATA  18 

REFERENCES  23 


APPENDIX  I  -  TECHNICAL  REPORT  "ANALYSIS  OF 
TIME-TO-EVENT  DATA  SUPPLIED 
BY  EGLIN  A.F.B. 

APPENDIX  II  -  INFERENCE  ON  THE  WEIBULL 
LOCATION  PARAMETER 

APPENDIX  III  -  INFERENCE  IN  WEIBULL  REGRESSION 


SKF  TECHNOLOGY  SEOv  CES 

SkF  \DbS'~:55  '  .C 

K.ING  OF  RUSSIA. 


AT80D048 


1.0  SINGLE  CLASSIFICATION  EXPERIMENTS 

A  paper  containing  a  comprehensive  treatment  of  sets  of 
two-parameter  Weibull  data  arising  from  single  factor  experi¬ 
ments,  was  published  in  the  Journal  of  Statistical  Planning 
and  Inference  (Vol.  3,  1979,  pp.  39  -  68).  A  second  paper, 
aimed  primarily  at  an  engineering  audience  was  prepared  and 
submitted  to  the  American  Society  For  Testing  and  Materials  (ASTM) 
for  presentation  at-  their  International  Symposium  on  Contact 
Rolling  Fatigue  Testing  of  Bearing  Steels  to  be  held  in  May  1981. 
The  paper  will  appear  in  the  Proceedings  of  that  meeting,  which 
is  planned  to  appear  as  an  ASTM  special  technical  publication. 

This  paper  extends  the  range  of  the  tabular  values  needed 
in  conducting  the  analysis  from  a  maximum  sample  size  of  n  =  10 
to  a  maximum  of  n  *  30.  Specifically,  new  tables  are  included 
for  n  *  15(5)30  with  censoring  amounts  rang'ing  from  r  *  5(5)n 
and  with  k  »  2(1)5  samples  per  experiment.  In  the  paper,  two 
illustrative  numerical  examples  are  used  to  illustrate  each  of 
the  analysis  procedures  viz.: 

(1)  Testing  the  equality  of  shape  parameters 

(2)  Testing  the  equality  of  scale  parameters 

(3)  Setting  confidence  limits  on  the  shape  parameter 

(4)  Setting  confidence  limits  on  the  tenth  percentile 

of  each  population  sampled,  and 

(5)  Conducting  a  single  range  multiple  comparison 
test  to  divide  the  population  into  groups. 
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The  data  used  for  illustration  were  the  results  of  roll¬ 
ing  contact  fatigue  tests  conducted  with  several  types  of 
steel  using  two  distinct  types  of  tester.  The  data  were  collect¬ 
ed  under  an  Air  Force  sponsored  program  with  one  set  of  tests 
conducted  at  Wright-Patterson  Air  Force  Base  and  the  other  at 
the  Pratt  and  Whitney  Division  of  United  Technology,  Inc. 

The  analysis  shows  that  the  life  ranking  of  the  steels 
■when  analyzed  as  a  set  was  the  same  for  both  test  devices. 

This  methodology  is  now  routinely  applied  within  SKF  to  analyze 
data  taken  for  industrial  and  DOD  sponsors  as  well  as  for  cor¬ 
porate  use. 

An  application  of  the  analysis  to  a  study  of  the  compara¬ 
tive  effect  of  six  bearing  greases  on  the  life  of  automotive 
wheel  bearings,  was  conducted  under  U.  S.  Army  Contract  No. 
DAAK70-77 -C-0034  and  appears  in  the  final  report  -  submitted  to 
the  U.  S.  Army  Mobility  Equipment  Research  and  Development 
Command  at  Fort  Belvoir,  Virginia.  The  report  is  SKF  Number 
AL78T022  entitled  "Performance  of  Automotive  Wheel  Bearing 
Greases . " 

2.0  LOCATION  PARAMETER  ESTIMATION  AND  INFERENCE 

During  this  contract  year,  work  on  inference  for  the  Weibull 
location  parameter  was  extended,  applied,  and  reported  in  full. 

Extensions  comprised  (1)  the  computation  of  additional 
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critical  values  of  the  statistic  used  for  testing  whether  the 
Weibull  location  parameter  exceeds  zero.  Values  are  now  avail¬ 
able  for  sample  sizes  ranging  from  10  to  40.  (2)  Additional 

computations  were  performed  for  the  purpose  of  comparing  the 
power  of  this  inferential  technique  to  the  Mann-Fertig  method. 
The  power  was  found  to  be  almost  identical.  (3)  An  analytical 
proof  was  developed  of  the  decreasing  monotonicity  with  X  of  the 

A  A  A 

ratio  w(X)  ■  6(ri)/8(r2)>  where  r^ <r2  and  where  g(r)  is  the  ML 
estimator  of  the  Weibull  shape  parameter  based  on  the  first  and 
r-th  order  statistics  when  an  amount  X  is  subtracted  from  each 
order  statistic. 


The  more  general  result  regarding  the  monotonicity  of 

«<* 

w(X)  when  f}(r)  is  defined  as  the  ML  estimator  based  on  all 
order  statistics  has  not  been  proven  but  has  repeatedly  been 
demonstrated  with  data  samples. 


I 

I 


-  Two  applications  of  the  methodology  were  made  using  actual 
data  samples  taken  from  two  diverse  areas  of  activity.  The 
first  application  was  to  a  data  sample  that  represented  time  to 
ignition  events  for  fuzing  devices.  The  possible  application  of 
the  three  parameter  Weibull  model  in  this  context  had  been  sug¬ 
gested  in  private  communication  with  Dr.  B.  Kurkji  an  when  he  was 
chief  mathematician  of  the  U.  S.  Army  Materiel  Command.  The 
model  is  attractive  for  this  application  because  if  it  fits,  the 


>*<“ 


v 


3 
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location  parameter  will  correspond  to  the  "safe"  time  prior 
to  which  an  armed  device  cannot  undergo  ignition.  Numerous 
contacts  were  made  at  Wright  -  Patterson  (W.  Romans  Air  Force 
Logistics),  at  Kirtland  (Neal  Chamblee) ,  and  at  Eglin  Air  Force 
Bases  to  discover  sources  of  time-to-event  data.  After 
discussions  with  staff  members  of  the  engineering  and  quality 
assurance  departments,  Mr.  Jasper  Glover,  head  of  the  Reliabili¬ 
ty  Department  at  Eglin,  finally  referred  us  to  Mr.  Charles  Yates 
of  his  department,  who  kindly  supplied  a  sample  of  29  time-to- 
event  data  unidentified  as  to  source.  The  data  were  analyzed 
and  showed  that  a  mixture  model  was  a  better  fit  to  the  data 
than  the  three  parameter  Weibull  model.  The  results  of  the 
analysis  were  presented  in  a  Technical  Report,  "Analysis  of 
Time-To-Event  Data  Supplied  by  Eglin  Air  Force  Base."  A  copy 
is  included  herein  as  Appendix  I.  Copies  of  the  report  were 
sent  to  Mr.  Yates  at  Eglin  and  on  his  recommendation,  to  Mr. 

T.  Mitchell  who  is  responsible  for  the  setting  of  safety  re¬ 
quirements.  The  results  were  also  discussed  with,  and  a  copy 
of  the  report  sent  to  a  Mr.  L.  Cox  of  the  Army's  Harry  Diamond 
Laboratories.  Mr.  Cox  is  concerned  with  the  performance  of 
fuzing  devices  developed  by  the  Army. 

The  Eglin  analysis  taught  that  the  methodology  is  useful 
for  distinguishing  two  superficially  similar  looking  types  of 
data,  namely:  samples  from  a  mixture  of  a  pair  of  2 -parameter 

4 
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Weibull  population.  To  make  this  distinction,  one  censors  a 
proportion  of  the  upper  end  of  the  data  and  reperforms  the  hy¬ 
pothesis  test.  If  a  significant  result  is  no  longer  so  after 
the  censoring,  the  mixture  model  is  assumed  to  obtain. 

A  second  application  was  made  using  data  taken  from  the 
results  of  a  foot  race.  The  object  was  to  determine  whether 
a  bound  on  human  performance  potential  could  be  found.  Again, 
it  was  found  that  a  Weibull  mixture  model  was  in  better  accord 
with  the  data  than  the  three  parameter  Weibull  model. 

The  extended  tabular  data  and  two  illustrative  examples 
discussed  above  have  been  described  in  a  paper,  "Inference 
on  the  Weibull  Location  Parameter,"  submitted  for  publication 
to  Technometrics.  A  verbal  presentation  having  the  same  title 
will  be  made  at  the  Joint  Meetings  of  the  American  Statistical 
Association  to  be  held  in  Houston,  Texas  in  August  1980.  A 
handout  synopsizing  the  talk  has  been  prepared  and  is  included 
herein  as  Appendix  II. 

3 . 0  WEIBULL  REGRESSION  OR  ACCELERATED  TEST  ANALYSIS 

During  a  previous  contract  year,  we  developed  methodology 
for  drawing  exact  inferences  in  this  setting:  (1)  Type  II  cen¬ 
sored  life  tests  are  conducted  at  various  levels  of  a  factor 
referred  to  as  a  stress.  (2)  At  each  stress  level,  the  life 
follows  a  two  parameter  Weibull  distribution  with  a  shape  para- 
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meter  8  that  is  invariant  with  stress  and  a  scale  parameter 
that  varies  inversely  with  a  power  y  of  the  stress.  Using  the 
distribution  of  certain  pivotal  functions  determined  by  Monte 
Carlo  sampling,  it  was  found  possible  to  set  confidence  limits 
on : 

(1)  The  exponent  in  the  relation  between  scale 
parameter  and  stress, 

(2)  The  Weibull  shape  parameter, 

(3)  A  percentile  of  the  life  distribution  at  any 
specific  stress. 

A  verbal  presentation  of  this  material  was  made  at  the  ASA 
joint  conferences  in  August  1979.  A  copy  of  the  handout  mater¬ 
ial  distributed  at  the  conference  is  given  in  Appendix  III. 

A  paper  was  prepared  describing  this  work  and  illustrating 
the  methodology  on  four  rolling  contact  fatigue  test  samples 
conducted  at  four  stress  levels.  Editorial  changes  to  the 
paper  were  made  in  the  current  contract  year  and  the  paper, 
entitled  "Confidence  Limits  for  Weibull  Regression  with  Censored 
Data,"  appeared  in  tie  IEEE  Transactions  on  Reliability,  Vol. 
R-29 ,  No.  2,  June  1980, 

A  question  raised  by  a  referee  regarding  the  goodness  of 
the  power  function  model  for  stress-life,  prompted  the  recogni- 

/s 

tion  that  one  could  perform  such  a  test  using  the  ratio  of  8(1) , 
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the  ML  shape  parameter  estimate  unconstrained  by  any  relation 
among  the  scale  parameters,  to  the  estimate  3  under  the  power 
function  constraint.  Accordingly  appropriate  code  was  added 
to  the  simulation  program  REGSIM  to  calculate  3(1). 

The  program  REGSIM  as  originally  configured,  calculates 
the  distribution  of  five  random  variables;  these  variables 

A 

being  specified  functions  of  the  ML  estimates  3,  and  the  p-th 

A 

quantile  Xp(S)  estimated  at  stress  S. 

The  program  REGSIM  has  now  been  modified  to  accommodate 
the  calculation  of  five  additional  random  variables  giving  a 
total  of  ten.  The  first  three  are  as  follows: 

(1)  3  /  3 

(2)  S( 1 ) / 6 

(3)  ( Y-Y)  .  S 

The  next  "k"  are  the  values  of  the  random  variable 

31n(xp/xp)  computed  at  the  k  stresses  ,  at"  which  life  tests  are  per- 
formed.  The  next  (7-k)  random  variables  are  the  values  of  3  . 

A 

In  (Xp/Xp)  at  (7-k),  other  specified  stresses. 

A 

The  distribution  of  the  first  random  variable  above,  3/3, 
is  needed  for  setting  confidence  limits  on  3.  The  second  is 
used  for  testing  the  adequacy  of  the  power  function  fit.  The 
distribution  of  the  third  random  variable  is  needed  for  calcu¬ 
lating  confidence  intervals  on  the  stress-life  exponent.  Finally 
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the  distribution  of  81n(kp/xp)  is  used  for  setting  confidence 
intervals  on  x  at  any  given  stress  level. 

Computer  runs  have  been  made  for  k  =  2,  with  life  tests 
at  stresses  sj  =  1.0  and  $2  =  1.2.  Supplementary  stresses 
were  taken  at  0.5  (0.1)  0.9.  For  k  =  3,  the  life  tests  were 
presumed  to  be  run  at  S  =  1.0,  1.1,  and  1.2  with  supplementary 
stresses  of  0.6(0. 1)0. 9. 

The  percentile  p  =  0.10  was  used  throughout  and  the  sam¬ 
ple  sizes  used  for  each  life  test  were  n  *  r  *  5,  10,  15,  20. 
Inasmuch  as  these  distributions  are  invariant  with  respect 
to  the  scale  of  the  stress  variable,  they  apply  when  the  stresses, 
in  whatever  physical  units  they  are  expressed,  are  proportional 
to  the  values  used  in  the  simulation  runs. 

A  short  paper  will  be  prepared,  aimed  at  a  user  audience, 
presenting  the  tables  and  illustrating  their  use. 

4.0  TWO-WAY  CLASSIFICATION  EXPERIMENTS 

The  likelihood  equations  for  a  general  two-way  factorial 
analysis  with  Weibull  response  have  been  derived.  There  are 
presumed  to  be  "a"  rows  and  "b"  columns,  representing  the  levels 
of  factors  A  and  B  respectively.  Five  separate  hypotheses  have 
been  considered  for  the  scale  parameter  applicable  when 

sampling  row  "i"  and  column  "j": 


'w'  —  w 
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Hi: 

nij 

*  ai  hj. 

H2: 

niJ 

“  ai  hi 

H3  * 

ni  j 

*  ai  n 

H4: 

nij 

=  bj  n 

H5: 

nij 

=  n 

a^,  bj  and  Cj^  are  multiplicative  row,  column  and  interaction 
effects  subject  to  the  constraints 

a  b  a  b 

"ai  =  tt  b j  =  "cij  =  it  cij  =  1 

i=l  j=2  i=l  j=2 


n  is  a  constant  "base  level"  scale  parameter  value.  is  the 

least  restrictive  hypothesis.  Under  each  cell  of  the  data 
layout  has  its  own  unique  scale  parameter  value.  is  the 

most  restrictive  hypothesis  under  which  all  cells  are  presumed 
ter  have  the  same  scale  parameter.  Under  H2  there  is  a  row 
and  column  effect,  but  no  interaction.  Under  H3  there  is  only 
a  row,  and  under  H* ,  only  a  column  effect. 

The  estimates  of  3  and  n x j  which  maximize  the  likelihood 
function  under  each  hypothesis  subject  to  the  row  and  column 
constraints  are  listed  in  Table  1.  They  represent  the  case 
where  n  items  are  tested  in  each  cell  until  the  first  r  fail. 
We  define  as  the  k-th  ordered  life  within  cell  (i,j). 
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TABLE  1.  ML  ESTIMATION  EQUATIONS  FOR  FACTORIAL 
EXPERIMENTS  UNDER  VARIOUS  HYPOTHESES 
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Define : 


Sij  =  2  In  [Xij(jc)  ] 
k=l 
a  b 

S-  *  *  *  SiJ 

J=l  1=1 
n  g 

Tij  =  2  *ij(k)ln  xij(k) 

ksl 

n  0 

vij  =  2  xij(k) 

ksl 


Vi.  =  (  TT  Vij)J 

j=l 

v.j  =  (  ?  Vij)- 

1=1 

V  s  (  v  n  ,r  \l/ab 


,!/b 

,1/a 


>1  ill  VU>' 


Likelihood  ratio  tests  can  now  be  constructed  to  test  the  more 
restrictive  of  these  hypotheses  against  less  restrictive  alter¬ 
natives.  For  example,  to  test  Hs,  the  hypothesis  that  is 
constant  over  i  and  j,  against  the  alternative  that  all 
differ,  one  would  calculate 


In  X  -  In  L(H5)  -  In  L(HL) 


where  LCH^)  denotes  the  likelihood  function  evaluated  using 
rj i j  and  B  estimated  by  the  methods  appropriate  for  hypothesis 
Hk. 

A  useful  sequence  begins  by  testing  Hg  against  Hi  to  assess 
the  hypothesis  Ho:Cjj  =  1,  i.e.  no  interaction.  If  H0  is  reject- 
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ed,  i.e.,  there  is  interaction,  no  other  tests  are  performed. 

If  there  is  no  interaction,  one  would  then  test  H-j  against 
H2  and  H4  against  H2  to  determine  respectively  whether  only 
column  or  row  effects  are  real. 

An  alternate  sequence  would  be  to  test  Hg  against  .  In 
this  case,  if  Hq  is  accepted,  no  further  tests  are  performed 
and  it  is  concluded  that  neither  row,  column  nor  interaction 
effects  are  significant.  If  Hq  is  rejected,  Hg  could  then 
be  tested  against  H2.  If  this  is  not  significant,  there  is 
interaction  and  testing  ceases.  If  it  is  significant,  there 
is  a  row  or  column  effect,  or  both.  One  then  tests  Hg  against 
H3  and  H4. 

In  either  sequence  of  tests,  H2  is  crucial  for  the  test 
of  interaction.  The  equation  for  estimating  the  shape  para¬ 
meter  under  H2  is  characteristically  different  from  the  esti¬ 
mating  equations  under  the  other  four  hypotheses  and  special 

A 

numerical  methods  will  need  to  be  developed  for  calculating  $2 • 

Under  H^,  the  estimation  equation  is  the  same  as  for  a 
single  factor  experiment  in  which  k  =  ab  tests  are  run.  Simi¬ 
larly,  under  Hg  the  estimating  equation  is  identical  to  that 
which  applies  to  a  single  sample  of  size  N  *  abn.  That  is,  one 
has  only  to  combine  the  data  in  all  ab  cells  of  the  design  into 
a  single  sample  and  estimate  the  shape  parameter  of  that  single 
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sample  to  obtain  85. 

The  estimating  equation  for  S3  is  of  the  same  form  as  the 
single  factor  experiment  with  k  »  a  and  with  the  data  in  each 

A 

row  combined,  i.e.  ignoring  the  columns.  Similarly  for  ^  the 
rows  are  ignored  and  the  data  within  each  of  the  columns  is  treated 
as  a  single  group  in  a  multiple  group  sample  with  k  *  b. 

Thus,  in  the  absence  of  interaction,  the  shape  parameter 
estimates  required  for  testing  row  and  column  effects  can  be  ob¬ 
tained  by  arranging  the  data  in  various  ways  using  only  the  soft¬ 
ware  for  ML  estimation  of  the  Weibull  shape  parameter  in  k  groups. 
The  appropriate  values  of  k  are  1,  a,  b,  and  ab.  Moreover,  with 
little  if  any  loss  of  power,  the  testing  can  be  based  on  just 
these  shape  parameter  estimates  to  avoid  the  need  for  additional 
software  to  calculate  the  likelihood  function. 

We  have,  accordingly,  generated  the  required  tables  for 
2X2,  2X3,  and  3X3  factorial  arrangements  with  the  sample 
sizes  n  and  censoring  number  r  tabled  below: 


Rows  X  Columns  n  r 


2X2  3  3 

4  4 

5  3 

5  5 

10  5 

10  10 
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Rows  X  Columns  n  r 


2X3  3  3 

4  2 

4  4 

5  2 

5  3 

S  5 

7  7 

3X3  3  3 

4  4 

5  5 


Under  the  assumption  that  interaction  is  negligible  (c^j  ■  l) , 
the  analysis  proceeds  as  follows: 

A 

1.  Combine  all  data,  calculate  85. 

2.  Treat  each  cell  as  a  separate  sample,  cal- 

A 

culate  81  • 

3.  Treat  each  row  as  a  separate  sample  (ignore 

A 

columns),  calculate  83. 

4.  Treat  each  column  as  a  separate  sample  (ig- 

A 

-  nore  rows),  calculate  84. 

A 

5.  Calculate  81/85.  If  greater  than  its  cri¬ 
tical  value,  row  effects  or  column  effects 
or  both  are  significant. 

A  A 

6.  Calculate  8^/83.  If  significant,  column 
effect  is  real. 

A  A 

7.  Calculate  S^/84.  If  significant,  row  effect 
is  real. 


XOGY 
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Table  2  shows  the  results  of  the  analysis  of  a  portion  of 
a  randomized  block  design  for  rolling  contact  endurance  test¬ 
ing  performed  by  Ku  et  al.  (1).  In  this  experiment,  it  was 
desired  to  determine  whether  there  was  a  difference  between 
two  oils  meeting  the  specifications  MIL-L-7808  and  MIL-L-23699 
with  respect  to  their  influence  on  fatigue  life  in  rolling 
contact.  Ten  specimens  were  run  to  failure  with  each  lubricant 
on  each  of  ten  test  machines.  We  have  arbitrarily  selected 
test  machines  Nos.  1  §  2  to  form  a  2  X  2  layout. 

The  values  shown  in  each  cell  are  (1)  the  ML  estimate  of 
the  tenth  percentile  x0.10  obtained  using  the  ten  data  values 
taken  at  the  conditions  corresponding  to  that  cell.  (2)  The  ML 
estimate  of  xo.lO  under  the  assumption  that  all  cells  have  a 
common  shape  parameter  and  (3)  the  ML  shape  parameter  estimate 
based  on  cell  data. 

A  A 

The  values  of  3^  and  65  are  shown  in  the  center  of  the 

A  A 

layout  in  Table  2.  83  and  84  are  shown  between  the  rows  and 

columns  respectively.  For  reference,  the  shape  parameter  es¬ 
timate  using  the  data  for  each  tester  are  shown  at  the  right 
hand  side  of  each  row.  The  shape  parameter  estimates  using  the 
combined  data  for  each  oil  are  given  below  each  column. 

To  test  the  homogeneity  of  shape  parameters,  an  assumption 
of  the  analysis,  one  computes  the  ratio  of  the  largest  to 
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TABLE  2.  ANALYSIS  OF  ROLLING  CONTACT  FATIGUE  DATA 


AT80D048 


smallest  of  the  cell  shape  parameter  estimates. 

w  »>  6.86/2.94  =  2.33 

The  critical  value  for  a  101  level  test  based  on  four 
samples  of  size  n  =  10  having  r  =  10  failures  is  (2). 

w0.90(n  =  10,  r  -  10,  k  -  4)  -  2.47 

The  hypothesis  of  homogeneous  shape  parameters  is  accept¬ 
ed  but  just  barely. 

Proceding  formally  we  form 

Bl/is  *  1-506 

This  is  substantially  greater  than  the  101  critical  value 

/V  /V 

(61/85)  q.90  =>  1.085  so  that  row,  column  or  both  effects  are 

significant.  Of  course,  an  interaction  effect,  suggested  by 
the  common  8  estimates  ofXQ.iQ,  also  would  contribute  to  a 
high  value  of  the  test  statistic. 

A  A 

To  test  the  difference  in  oils,  we  calculate  81/83  *  1.083. 

A  A 

The  10%  level  critical  value  is  (81/83)  0.90  =  1*055.  Thus,  the 
lubricant  difference  is  virtually  significant. 

To  assess  the  difference  between  testers,  we  compute 

A  A 

81/84  ■  1.47.  This  greatly  exceeds  the  critical  value  1.085, 
suggesting  a  strong  tester  effect.  This  effect  is  not  suggested 
by  the  common  8  shape  parameter  estimates  of  xq.10  These  es- 
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tiraates  in  fact  suggest  an  interaction  effect. 

The  marginal  shape  parameter  estimates  suggest  an  inter¬ 
action  as  well,  for  if  the  lube  effect  were  the  same  with  each 
tester,  the  shape  parameter  estimates  within  each  tester  would 
be  homogeneous.  The  ratio  of  the  tester  shape  parameter  esti¬ 
mates  is  5.54/3.48  3  1.59.  This  exceeds  the  critical  value 
w0.9o(20,20,2)  for  a  10%  level  test  of  the  homogeneity  of  two 
shape  parameters  based  on  censored  samples  of  size  20. 

This  effect  is  also  ascribable  to  the  low  shape  parameter 
estimate  for  the  Tester  No.  2  -  MIL-L-7808  cell. 

We  conclude  that,  as  in  the  analysis  of  variance  for 
normal  distribution  theory,  inhomogeneous  variance  (shape 
parameter)  can  cause  specious  results.  A  second  application 
was  made  to  rolling  contact  fatigue  data  taken  at  two  loads 
and  with  two  radii  of  curvature  (3).  These  data  also  exhibited 
a  nonhomogeneous  shape  parameter. 

A  paper  will  be  prepared  on  the  analysis  of  randomized 
block  designs  with  Weibull  response. 

5.0  GROUPED  DATA 

Computer  program  WEIBSIM  for  simulating  sets  of  Weibull 
distributed  data  has  been  modified  to  form  a  new  program  GROUPSIM 


r<~ 


\  n  ,  r*> - - ,  «  /- 

<  <V  v. 
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for  the  analysis  of  grouped  Weibull  data. 

The  program  generates  one  or  more  ungrouped  samples,  allo¬ 
cates  each  failed  item  to  an  interval,  replaces  its  actual  value 
by  the  cell  midpoint  value,  and  then  performs  conventional  max¬ 
imum  likelihood  estimation  using  the  data  thus  modified. 

The  program  assumes  type  II  censoring  at  the  r-th  failure.  The 
implication  is  that  testing  stops  when  a  pre-established  failure 
occurs.  In  actual  testing,  additional  failures  could  occur 
prior  to  the  end  of  the  interval  containing  the  r-th  failure. 

The  intervals  have  been  chosen  logarithmically.  An 
initial  interval  'DELT'  is  input  to  the  program  along  with  a 
factor,  'FAC', 

The  first  interval  extends  from  0  to  DELT.  The  second 
interval  extends  from  DELT  to  DELT  X  FAC.  The  terminus  of  the 
1-th  interval  is  calculated  as 

A (I)  =  DELT  X  (FAC)1'1 

or  recursively  as 

A  (I)  =  FAC  X  A  (1-1) 

The  program  samples  from  a  Weibull  population  having  shape 
parameter  8*1  and  a  p-th  percentile  Xp  =  1.0  for  specified  p. 

The  simulation  results  apply  to  any  2  parameter  Weibull  distri- 


19 


AT80D048 


bution  if  the  terminus  x^  of  the  i-th  interval  satisfies 

(xi/xp)3  =  A  (I) 

For  the  distribution  used  in  the  simulation,  the  1st  and 
99-th  percentiles  are  0.0954  and  43.70  when  using  p  =  0.10. 

GROUPSIM  was  used  to  determine  the  distribution  of  the 

/V  *  ^ 

pivotal  quantities  S/ S  and  Sin  (xq . io/*o . 10)  for  a  single  un¬ 
censored  sample  of  size  n  =  30,  using  DELT  =  0.1  and  FAC  =  2.0, 
1 . S ,  and  1.2. 

The  5-th,  50-th  and  95-th  percentiles  are  tabled  below  along 
with  the  corresponding  ungrouped  values. 


B/8 

§ln(x0 

h-» 

O 

X 

o 

o 

v — * 

0.05 

0.50 

0.95 

0.05 

0.50 

0.95 

FAC  *  2.0 

0.900 

1.056 

1.269 

-0.346 

0.138 

0.742 

FAC  *  1.5 

0.927 

1.073 

1.267 

-0.306 

0.132 

0.672 

FAC  *  1.2 

0.951 

1.080 

1.239 

-0.234 

0.131 

0.557 

UNGROUPED 

0.826 

1.057 

1.335 

-0.567 

0.0536 

0.915 

To  assess 

the  sample 

size  effect  the 

values  for 

n  =  r  =  5 

and 

n  =  r  ■  50 

with  DELT 

=  0.1  and 

FAC  = 

1.5  are  tabled  below. 
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»* 


S/8 

8lnCx0>1o/x0>10) 

0.05 

0.50 

0.95 

0.05  0.50  0.95 

n  =  r  *  5 

0.768 

1.066 

1.654 

-0.869  -0.129  0.913 

UNGROUPED 

0.680 

1.235 

2.815 

-1.142  0.447  4.445 

n  *  r  *  50 

0.937 

1.060 

1.213 

-0.252  0.121  0.564 

UNGROUPED 

0.852 

1.018 

1.235 

The  following  effects  are  observed: 

(1)  There  is  consistently  less  variability  in  the  grouped 
data  results  than  the  associated  ungrouped  values. 

(2)  For  fixed  sample  size,  the  variability  decreases  with 
the  interval  width  as  expected,  but  does  not  appear 
to  be  converging  toward  the  ungrouped  results. 

(3)  The  difference  between  grouped  and  ungrouped  percen¬ 
tage  decreases  with  sample  size,  i.e.  convergence  with 
sample  size  appears  to  take  place. 


I 

I 

I 


Superficially  it  appears  that  grouping  the  data  results  in  great 
er  precision  in  estimating  the  parameters  than  ungrouped  data, 
a  counter  intuitive  result.  The  grouping,  however,  assumes  per¬ 
fect  information  regarding  the  parameters  for  the  purpose  of 
standardizing  the  intervals.  That  is,  the  values  of 
(x^/xq.iq)®  defining  the  interval  end  points  is  assumed 
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known,  whereas,  in  fact,  only  the  are  known.  It  is  recommend¬ 
ed  that  this  effect  be  examined  in  future  studies.  It  is  further 
recommended  that  the  results  given  herein,  obtained  by  conventional 
ML  estimation  using  adjusted  data,  be  compared  in  future  studies 
to  those  obtained  with  direct  ML  estimation  using  a  grouped  data 


formulation. 
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ANALYSIS  OF  TIME-TO-EVENT  DATA 
SUPPLIED  BY  EGLIN  AFB 


1 . 0  Introduction  and  Summary 

A  methodology  has  recently  been  developed  under  the  sponsor¬ 
ship  of  the  Air  Force  Office  of  Scientific  Research  for  the 
unbiased  point  and  interval  estimation  of  the  location  parameter 
of  a  three-parameter  Weibull  distribution. 

It  has  been  suggested  [1]  that  time-to-event  data  of  the 
type  encountered  in,  fuzing  mechanisms  may  follow  a  three-parameter 
Weibull  distribution.  If  this  is  so  the  location  or  threshold 
parameter  of  the  distribution  represents  the  "safe"  time  prior 
to  which  the  ignition  event  can  not  occur. 

In  this  context  a  lower  confidence  limit  on  the  location 
parameter  represents  a  quantifiably  conservative  estimate  of  the 
"safe  time"  for  the  device. 

This  report  describes  the  analysis  of' a  sample  of  29  uniden¬ 
tified  time-to-event  observations  supplied  by  Eglin  Air  Force  Base. 

Under  the  assumption  that  the  data  are  drawn  from  a  three- 
parameter  Weibull  population,  a  median  unbiased  estimate  and  a 
lower  95%  confidence  limit  for  the  location  parameter  have  been 
calculated  using  computer  program  "LOCEST"  implementing  the  metho¬ 
dology  referred  to  above.  Subtracting  the  median  unbiased  location 
parameter  estimate  from  each  event  time  and  regarding  the  data  thus 
adjusted  as  a  two-parameter  Weibull  sample,  computer  program 
"WEIBEST"  was  used  to  estimate  the  shape  and  scale  parameter  by 
the  method  of  maximum  likelihood.  A  chi-square  goodness  of  fit 
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test  was  then  performed  and  indicated  that  the  fitted  three- 
parameter  Weibull  distribution  did  not  adequately  describe  the 
data.  Instead,  the  data  appear  to  consist  of  a  mixture  of  two 
Weibull  distributions  with  one  population,  representing  roughly 
86%  of  the  data,  having  a  shape  parameter  of  6  *  18.6  and  a  second 
population  representing  14%  of  the  observations  and  having  a  shape 
parameter  estimated  graphically  to  be  8  *  1.4. 

The  data  thus  support  the  assumption  that  the  population 
from  which  the  sample  was  drawn  consist  of  a  mixture  of  effective 
items  having  a  two-parameter  Weibull  distribution  and  a  10-20% 
subpopulation  of  "duds"  for  which  the  event  occurs  at  random 
intervals  in  accordance  with  a  poisson  process.  To  estimate  the 
"safe"  time  associated  with  some  arbitrary  low  event  probability, 
the  long  lived  items  should  be  censored  prior  to  estimating  the 
Weibull  parameters.  If  this  is  not  done  the  Weibull  shape  para¬ 
meter  for  the  effective  items  will  be  underestimated  and  overly 
conservative  safe  lives  will  be  computed. 

Section  2.0  of  this  report  describes  the  methodology  for 
location  parameter  estimation.  Section  3.0  gives  the  details 
of  the  analysis  of  the  Eglin  time-to-event  data. 
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2 . 0  Methodology 

2 . 1  The  Weibull  Distribution 

The  three  parameter  Weibull  distribution  has  the  cumulative 

form: 

Prob(X  <  x]  »  F(x)  *  1  -  exp  [-fx-yVnJ®]  ;  x  >  y  (1) 


where  y  *  location  parameter 
n  ■  scale  parameter 
8  *  shape  parameter 


The  two-parameter  Weibull  distribution  is  the  special  case 
of  Equation  (1)  in  which  the  location  parameter  y  -  0. 

2 . 2  Graph  of  Two  and  Three  Parameter  Weibull  Functions 

For  the  two-parameter  Weibull  distribution  it  is  readily 
shown  that 

y(x)  *  inln(l/(l-F(x)  )  *6  In  x  -  Q  In  r\  ( 2 ) 

Thus,  in  the  two-parameter  case  y(x)  is  a  linear  function  of 
in  x  having  slope  8  and  intercept  For  the  three-parameter 

case 

y(x)  =  in  in  (l/(l-?(x))  =  3  In(x-y)  -  3  lr.  n  ;  x  >  y  (!) 


The  slope  of  a  plot  of  y(x)  against  Inx  in  the  three-parameter 
case  is 


du(x)  m  8a 
dln(x)  x-y 


;  x  >  y 


(4) 
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The  slope  is  infinite  at  x  »  y  and  decreases  monotonically 
thereafter  to  an  asymptote  of  8.  Figure  1  is  a  sketch  of  y(x) 
plotted  against  Inx  for  y  =  0  and  y  >  0. 

2 . 3  Graphical  Estimation  of  Shape  Parameter 

Let  denote  the  ordered  observations  in  a 

l  u  O  71 

random  sample  of  size  n  drawn  from  a  two  or  three-parameter  Weibull 
distribution.  An  estimate  of  may  be  calculated  for  each 

of  the  ordered  observations  using  any  of  the  various  choices  of 
plotting  position.  A  common  choice  is: 

F(X‘)*i/(n  +  l)  (5)' 

u 

An  estimate  y(x.)  may  then  be  computed  by  substituting 

u 

a 

F(xi)  into  Equation  (2). 

If  the  sample  is  drawn  from  a  two-parameter  Weibull  distri- 
bution  y(x.)  will  tend  to  plot  against  Inx .  as  a  straight  line 
with  slope  8.  If  y >0,  i.e.  the  population  is  a  three-parameter 

A 

Weibull  distribution,  y(x.)  will  tend  to  be  a  concave  function 
of  inx^  approaching  a  constant  slope  8  for  large  x^  values. 

Figure  2  shows  how  a  plot  of  y(x^)  vs.  Inx ±  might  appear 
for  a  sample  drawn  from  a  three -parameter  Weibull  distribution. 

If  these  data  were  regarded  as  a  two -parameter  Weibull  sample 

A 

a  graphical  estimate  of  the  shape  parameter  6^  could  be  found  as 
the  slope  of  the  straight  line  that  best  fits  the  complete  data 
sample . 

If  only  a  subset  of  the  smallest  ordered  values  were  used 
in  graphically  estimating  the  shape  parameter,  the  estimate 

A 

3-  would  be  obtained.  For  three  parameter  Weibull  data  3^  will 
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tend  to  exceed  QA.  On  the  other  hand  when  the  sample  is  drawn 

A  A 

from  a  two-parameter  Weibull  population  (y =0),  8 £  and  84  will 
be  comparable. 

2 . 4  Maximum  Likelihood  Estimation  of  Shape  Parameter 

Rather  than  graphical  estimation  we  consider  maximum  likeli 
hood  (ML)  estimation  of  8. 

For  a  sample  of  size  n  censored  at  the  r-th  ordered  observa 
tion  xp,  the  ML  estimate  of  3  for  a  two-parameter  Weibull  distri 
bution  is  the  solution  of  the  nonlinear  equation: 

A  ji  p  A  A  p  A  A 

1/ 8  +  Z  logxy/r  -  (  Z  x^  logx^+x^(n-r)x  )/(  Z  x$  +  (n-r)x®)  =  0  (6) 

t»2  i-1  r  i-1 

As 

It  has  been  shown  (cf.  McCool  [2])  that  8/8  is  a  pivotal 
function,  i.e.  it  follows  a  distribution  that  depends  on  n  and 
r  but  not  on  the  underlying  Weibull  population  parameters. 

A 

Denoting  the  solution  of  Eq.  (6)  as  8^  and  the  solution  of 

/\  /S  /N 

Eq.  (6)  with  r2<r  as  S>L,  the  distribution  of  w  =  8^/8 £,  will 
depend  only  upon  n,  and  r  when  the  underlying  distribution 

is  indeed  of  the  two-parameter  Weibull  form.  When  the  underlyin 
distribution  is  the  three-parameter  Weibull  the  mean  value  of 

A 

8r  will  increase  proportionally  more  than  the  mean  value  of  3  i . 

u 

With  the  percentiles  of  w  determined  by  Monte  Carlo  samplin 
for  specified  r2  ,  r  and  n,  one  may  reject  the  hypothesis  that 
y  m  o  at  the  100a%  level  if 
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2 . 5  Interval  Estimation 

Given  that  a  random  variable  x  is  drawn  from  a  three-parameter 
Weibull  population  having  location  parameter  y  the  transformed 
variable  y  =  x-y  will  follow  a  two-parameter  Weibull  distribution 
with  the  same  scale  and  shape  parameters  as  the  three-parameter 
distribution.  Thus,  if  y  is  subtracted  from  the  observed  data 

A  A 

prior  to  calculating  84  and  8^  from  Eq.  (6),  the  resulting  ratio, 
denoted 

w(y)  =  ZA/SL  (8) 

will  follow  the  null  distribution  of  w  determined  by  Monte  Carlo 
sampling  from  a  two-parameter  Weibull  population  for  given  values 
of  n,  Vj  and  r. 

We  may  thus  write  the  100(l-a)%  probability  statement 
Prob  [w(y)  <  Uj-a.]  ~  1  ~  a 

We  also  need  the  fact,  heuristically  suggested  by  the  analogy 
to  graphical  estimation,  that  if  an  amount  X  is  subtracted  from 
each  observation  in  a  given  sample  prior  to  calculating 

A  A  A  A 

84  and  8^ ,  w(\)  -  8^/3 £  will  be  a  decreasing  function  of  X. 
Accordingly,  we  may  invert  the  inequality  of  Eq .  (9)  to  give  a 
100(l-a)%  lower  confidence  limit  for  y,  i.e. 


y  >  w~ 1  (10) 
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3.0  Analysis  of  Eglin  Data 


3.1  Analysis  as  a  Two-Parameter  Weibull  Sample 


Table  1  is  the  output  of  computer  program  W'EIBEST  which  was 
applied  to  the  raw  Eglin  data.  The  tabular  output  at  the  top  of 
Table  1  is  a  sorted  list  of  the  29  observations  of  time-to-event 
in  seconds. 

The  fi^st  line  below  the  sorted  times  gives  the  maximum 
likelihood  estimates  of  the  tenth  and  50-th  percentiles  (desig¬ 
nated  L-j^q  and  L.5Q,  respectively)  and  the  Weibull  shape  parameter 
6,  computed  under  the  assumption  that  the  data  were  drawn  from 
a  two-parameter  Weibull  distribution.  Subsequent  lines  in  Table 
1  give  lower  and  upper  90%  confidence  limits  and  T"cdi'!n  unbiased 
estimates  of  L10,  L5Q  and  g. 

Figure  3  shows  a  probability  plot  of  the  data  using  scales 
on  which  two  parameter  Weibull  samples  tend  to  plot  as  a  straight 
line.  The  fitted  two-parameter  population  is  shown  as  a  solid 
straight  line  and  is  clearly  a  poor  fit  to  the  data.  The  two 
dotted  straight  line  segments  fitted  to  the  data  are  discussed 
further  below. 

3.2  Analvsis  as  a  Three-Parameter  Weibull  Sample 


Figure  4  is  a  plot  of  the  function  u ( X )  computed  from  Eq. 

(8)  for  positive  X  values.  In  calculating  this  plot  was  taken 
as  5  and  r2  as  29.  The  plot  decreases  with  X,  approaching  a  ver¬ 
tical  asymptote  as  X  approaches  the  first  order  statistic 


TABLE  1 
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WEIBEST  OUTPUT  RAW  DATA 


TIME-TO-EVENT  DATA  FROM  EGLIN  AFB 


Group  No.  1  Lives 


5.2200 

6.0700 

6.5500 

5.4100 

6.1500 

6.5600 

5.7300 

6.2600 

6.5900 

5-7800 

6.3100 

6.6000 

5-7900 

6.3400 

6.7300 

5.8700 

6.3500 

7.1300 

5-9200 

6.4300 

7.8300 

5.9800 

6.4400 

9.2800 

6.0000 

6.4500 

10.2800 

6.0300 

6.4600 

ilfl. 

L^0_ 

BETA 

0.4556E  01 

0.6501E  01 

0.5299E  01 

LCL  L^q 

MED  L10 

UCL  L10 

0.3748E  01 

0.4524E  01 

0.5113E  01 

LCL  L^q 

MED  L50 

UCL  L50 

0.6024E  01 

0.6499E  01 

0.6912E  01 

MED  BETA 


LCL  BETA 

0.3893E  01 


0.5125E  01 


UCL  BETA 
0.6592E  01 


TIME-TO-EVENT  DATA  FROM  EGL.IN  AFB 
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The  median  unbiased  estimate  yQ  ^is  shown  to  be  the  X  value 
associated  with  w(\)  -  w.  and  the  95%  lower  confidence  limit 

A 

Y o.OS  *s  the  *  value  corresponding  to  wQ  gg. 

Note  that  if  w(0)<  wQ  gs  ,  a  positive  value  of  y q,qs  cannot 
be  found. 

Computer  program  LOCEST  calculates  yQ  SQ  and  yQ  Q5  for  spe¬ 
cified  values  of  w .g5  and  w.  Sq  using  a  golden  section  search  tech¬ 
nique.  The  appropriate  values  found  from  Monte  Carlo  simulation 
for  n*30,  i,2s‘s*  r  a 90  are 

w0 . SO  *  1 ' 294 
^  9  5  ~  3.279 

The  values  for  n-29  ,  r=2P  are  not  likely  to  differ  substan¬ 
tially  from  these.  Table  2  shows  the  LOCEST  output. 

The  Weibull  shape  parameter  considering  24  items  censored 
at  the  5-th  smallest  time  is 
lL  -  27.9 

Using  all  the  data  the  estimate  is 
-  S.30 

A  A 

The  ratio  B^/8^  55  S.27  corresponds  to  w(o)  and  since  w(o)  exceeds 
w0.50  and  W0.9S>  positive  values  of  both  ygm5Q  and  y g.os  ma>r 
calculated. 

These  values  are  yg  =  3.88 
and  ^O.SO  a  s'09 

A 

The  two  parameter  estimates  of  ioaLio~'i ,  so  and  s0.50 

=*Lso-y.SO  are  0.364  and  1.243,  respectively.  The  estimated  shape 
parameter  from  the  adjusted  data  is  1.53. 
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LOCEST  OUTPUT 

TIME-TO-EVENT  DATA  FROM  EGLIN  AFB 
WEIBULL  LOCATION  PARAMETER  ESTIMATION 

SAMPLE  SIZE ,N=  29 

TRUNCATION  NUMBER  FOR  LOCATION  PARAMETER  ESTIMATION , Rls  5 
NUMBER  OF  FAILURES, R=  29 
W50=  1.291* 

W95=  3.279 

BETA  HAT'(R1)=  27-940 

BETA  HAT(R)  =  5.299 

BETA  HAT( R1 ) /BETA  HAT( R ) =  5.272 

MEDIAN  UNBIASED  ESTIMATE  OF  GAMMA=  5.091 

LOWER  95%  CONFIDENCE  LIMIT  FOR  GAMMA=  3-881 

ADJUSTED  MAXLIKE  ESTIMATES  OF  10-TH  AND  50-TH  PERCENTILE 

X0.10=  0.364 

XO . 50=  1.243  ADJUSTED  BETA  r  1.533 
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Figure  5  is  a  probability  plot  of  the  data  after  adjustment 

A 

by  subtraction  of  Y Q  s0,  with  the  fitted  distribution  shown  as 
a  solid  straight  line.  A  chi-square  goodness  of  fit  test  was 
applied  and  led  to  rejection  of  the  hypothesis  that  the  fitted 
distribution  actually  represented  the  data  (cf.  Appendix).  These 
dashed  line  segments  emphasize  that  the  behavior  in  the  two  tails 
is  inconsistent  with  the  fitted  distribution. 

Returning  to  Figure  3,  we  note  that  the  two  dashed  line  seg¬ 
ments  together  fit  the  observed  data  very  well  and  suggest  that 
the  data  may  be  a  mixture  of  two  Weibull  populations;  one  popula¬ 
tion  having  a  high  shape  parameter  value  and  low  mean  time-to-event 
and  a  second  population  having  a  much  lower  shape  parameter  value 
and  a  high  mean  time-to-event.  For  this  data  sample  86%  of  the 
sample  belongs  to  the  first  population. 

To  estimate  the  parameters  of  the  first  population,  the  data 
were  censored  at  the  24-th  event  time  and  WEIBEST  was  rerun.  The 
results  are  given  in  Table  3. 

The  shape  parameter  is  estimated  to  be  18.6.  This  is  much 
higher  than  the  value  5.3  shown  in  Table  2  based  on  all  the  data. 

Figure  6  is  a  probability  plot  of  the  censored  data  and  con¬ 
firms  the  good  fit  of  the  two-parameter  V.reibull  population  to  the 
bulk  of  the  early  events. 

A  graphical  estimate  of  the  shape  parameter  for  the  long 
event  time  items  is  1.4.  This  is  consistent  with  a  shape  para¬ 
meter  of  1.0  which  suggests  that  the  long  event  time  population 
may  have  an  exponentially  distributed  time  between  failures 

characteristic  of  a  poisson  process  governing  the  occurrence  of 
events.  This  suggests  that  the  late  events  correspond  to  a  sub- 
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TABLE  3 

WEIBEST  OUTPUT  FOR  CENSORED  DATA 


EGLIN  DATA  CENSORED  AT  24-TH  ORDERED  OBSERVATION 


Group  No.  1  Lives 


5.2200 

6.0700 

6.5500 

5.4100 

6.1500 

6.5600 

5.7300 

6.2600 

6.5900 

5.7800 

6.3100 

6.6OOOS 

5.7900 

6.3^00 

6.6OOOS 

5.8700 

6.3500 

6.6000S 

5.9200 

6.4300 

6.6000 

5.9800 

6.4400 

6.6000S 

6.0000 

6.4500 

6.6000S 

6.0300 

6.4600 

LiQ 

Jacl- 

BETA 

0.5695E  01 

0.6301E  01 

0.1863E  02 

LCL  Li o 

MED  L^o 

UCL  L10 

0.5366E  01 

0.5677E  01 

0.5870E  01 

LCL  L50 

UCL  Lqn 

0.6166E  01 

C.6303E  01 

0.6436E  01 

LCL  BETA 

MED  BETA 

UCL  BETA 

0.1296E  02 

0.1782E  02 

0.2377E  02 
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population  of  "duds"  wherein  the  event  is  triggered  by  some  sort 
of  random  shock  rather  than  by  the  design  mechanism. 

3 . 3  Discussion  of  Data  Analysis 

The  data  do  not  confirm  a  three-parameter  Weibull  model  im¬ 
plying  a  "safe"  time  prior  to  which  the  "event"  cannot  occur. 

A  "safe"  time  must  therefore  be  defined  as  the  time  associated 
with  some  arbitrary  low  probability  that  the  event  will  occur 
prior  to  it.  Because  the  data  suggest  a  mixture  of  Weibull  models 
the  direct  use  of  a  two-parameter  Weibull  model,  as  shown  in 
Figure  3,  will  result  in  overly  conservative  safe  time  estimates. 
This  may  be  overcome  by  censoring  15-20%  of  the  long  time  events 
which  correspond  to  a  population  of  "duds"  that  are  mixed  with 
the  effective  items. 
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APPENDIX 

GOODNESS  OF  FIT  TEST  FOR  THREE -PARAMETER  WEIBULL 

From  Tables  2  and  3  the  fitted  three-parameter  Weibull  CDF 
has  the  equation 

F(s)  -  1  -  exp  -  0.10526  2  ’  524  (A.l) 

where 

t  *  x  -  S .  0 9  (A. 2) 

We  divide  the  t  axis  into  S  intervals  each  having  a  20%  occur 
rence  probability  by  calculating  the  percentiles  tQ  2o> 
t0.40 »  *0.80  and  *.80  where  *p  satisfies 

p  ~  1  -  exp  -  .10526  [ vp/.2642 ]  (A. 3) 

The  expected  number  of  observations  in  the  i-th  interval  is 

a  .  »  n  x  p  *  29  x  .2(7  »  5.8  (i  *  1,...5)  (A.  4) 

The  following  Table  shows  for  each  interval  the  observed 

number  of  observations  in  the  interval,  o^,  the  expected  number, 

2 

e £ ,  and  the  square  of  the  differences  (o^-  e ,■). 


Interval  on  t 

Q; 

gy 

(  O  y  -  g  : 

0  -  0.5941 

2 

5.8 

14.44 

>  0.5941  -  1.020 

9 

5.8 

10.24 

>  1.020  -  1.492 

11 

5.8 

27.04 

>  1.492  -  2.154 

4 

5.8 

3.24 

>  2.154 


00 


3 


S .  8 


7.84 
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Under  the  hypothesis  that  the  data  were  drawn  from  the  fitted 
distribution,  the  quantity 

5  2 

u-  l  (oi  -  &i)  /ei  =  S6.2/S.8  *  9.89  (A.  5) 

i-1 

asymptotically  follows  a  x 2  distribution  with  m-p-1  degrees  of 
freedom  where: 

m  *  no.  of  intervals 

p  *  no.  of  parameters  estimated  by  the  method 
of  maximum  likelihood 

This  asymptotic  result  is  generally  believed  to  be  applicable 
if  *  S .  In  the  present  case  m=5  and  p=3  so  that  u  will  be 
approximately  distributed  as  X2(l)  under  the  null  hypothesis. 

p 

From  tables  of  the  x  distribution  we  find 

Since  u  *  9.69  >  3.84  the  null  hypothesis  is  rejected  at  the  5% 
significance  level. 
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INFERENCE  ON  THE  WEIBULL  LOCATION  PARAMETER 


J.  I.  McCool 
SKF  Industries,  Inc. 
King  of  Prussia,  Pa.  19406 


1.  HYPOTHESIS  TEST;  2  vs.  3  PARAMETER  WEIBULL 

The  two-parameter  Weibull  distribution  has  CDF 

8 

F(x)  =  l-exp-[x/n] 

The  maximum  likelihood  estimate  of  the  shape  parameter  B 
calculated  from  the  ordered  observations  xi  in  a  random  sample 
of  size  n  type  II  censored  at  the  r-th  failure  is  denoted  &(r) 
and  is  the  solution  of: 

^  A 

1/8  +  5  ln(xi)/r-  ?  Xi^ln(xj)/  ?  x*®  =  0 

iol  i=l  i= 1 

where  x^  =  xr,  i  >  r. 

Define  the  random  variable 

A 

w  =  8  ( r i )/  8(r) 

where  ri  <r.  The  distribution  of  w  depends  only  on  rj,  r 
and  n. 


When  sampling  from  the  three-parameter  Weibull  distri¬ 
bution  with  location  parameter  y  >  0,  w  becomes  stochastically 
larger.  The  acceptance  region  for  a  100a£  level  test  of  Hq :  y= 
0  against  against  Hi:  Y  >0  is: 


SKF 


■|\ 


Z  '< 


1 
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8  (r-| )/  S  (r)  <  wi_  a 

Percentage  points  found  by  Monte  Carlo  sampling  are 
given  in  Table  1.  In'  studies  with  6  =  1  the  choice  r-|=5  was 
found  to  be  nearly  optimum  for  all  n  and  r  in  the  range 
represented  in  Table  1. 


2.  INTERVAL  ESTIMATION 

Subtract  Y(>0)  from  each  ordered  observation  in  a 
random  sample  from  a  three-parameter  Weibull  distribution. 
Define 

w(X  )  =  8(r-] )/  8(r) . 

w(X )  is  a  decreasing  function  of  X  and  w(X=y)  follows  the 
same  distribution  as  w  calculated  from  two-parameter  Weibull 
samples. 


Inverting  the  statement:  Prob[w(y)<  wi_a]  =  l-ct  gives 
the  lower  100(l-a)$  confidence  limit: 

Y  >  w_1  (w-|_a)  • 


3. 


EXAMPLE 

The  times  to 
are  as  follows: 
5.22,  5.41,  5.73, 
6.07,  6.15,  6.26, 
6.55,  6.56,  6.59, 


an  ignition  event  f 

5.78  5.79,  5.87,  5 
6.31,  6.34,  6.35,  6 
6.60,  6.73,  7.13,  7 


29  fuzing  devices 


32, 

5.98, 

6.00, 

6.03, 

43, 

6.44, 

6.45, 

6.46, 

83, 

9.28, 

10.28 

Using  ri=5  and  r=29  gives 

A  A 

8  (rj)  .-  27-9  and  8(r)  =  5-30 
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A  A 

B(ri)/8(r)  =  5.27  >  3-279  =  110.95(^  =  5  ,r=30) 

The  two-parameter  Weibull  hypothesis  is  rejected  in 
favor  of  the  three-parameter  alternative. 

The  median  unbiased  estimate  of  the  location  parameter, 
y0.50  is  found  by  solving 

w(  Yo . 5o)  s  w0.50  =  1-294 

to  be 

Y0.50  =.  5.09 

A  lower  95i  confidence  limit  Y 0.50  was  found  by  solving 

A 

W(Y0.05)  =  w0.95  =  3-279 

to  be 

a 

Y0 . 05  s  3-88 
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THE  MODEL 

A  random  variable  x  is  presumed  to  follow  a  two-parameter 
Weibull  distribution  with  shape  parameter  8  and  scale  parameter 
n  that  varies  as  an  inverse  power  of  a  positive  deterministic 
external  variable  S,  generically  termed  a  stress.  That  is, 

F(x|S)  -  1  -  exp  -  [x/n(S)  ] 8  (1) 

n(s)  *  n0  S*Y  C2) 

with  n  ,  8,  y,x>  o 

As  a  consequence  of  Cl)  and  (2)  the  p-th  quantile  at  stress  S  is 
Xp(S)  -  (-in  (1-p) ) 1/6  •  n(S)  (3) 

MAXIMUM  LIKELIHOOD  CML)  ESTIMATION 

A  life  test  is  carried  out  at  each  of  k  stress  levels  denoted 
S1 ’  ^2 ‘ \ •  At  the  i-th  stress  level  ni  items  are  tested  until 
the  first  r^  ordered  failures  are  observed  (type  II  testing). 

The  ML  estimates  of  Y  and  8  are  calculated  as  the  simultaneous 
solution  of  the  following  equations. 


AL79P02  2 


AA  H  •  A 

ys  / 


Z  Z  logx.j---.  Z  Z  x^,...  logXjn-) 

1  ^  i-l  i*l  1UJ  i-l _ _  j-1  UJ  1U  J 


Z  r. 


z  sY0  rix*  , 

i-l  1  j-1  1(^ 


where  is  j“th  ordered  life  in  the  i-th  sample. 


.The  ML  estimate  of  n  is 

o 


V  ^4  y  A  If 

tr  i1  A OV)/:  n 

i-l  j-1  Si  Y  i*l 


}  1/6 


The  ML  estimates  of  nCs)  or  xpC.S)  are  given  by  Eqs.  (1)  and 
(2)  on  substituting  the  ML  estimates  of  y,  6  and  n. 


PIVOTAL  FUNCTIONS 

The  following  functions  are  "pivotal",  i.e.  they  follow 
distributions  that  do  not  depend  on  the  population  values  of 
the  parameters: 

q  -  6/  8  0 

A  A 

w*  -  (y  -  y)  6  CS 

U*  -  a  in  lipCS)/xpCS)]  C£ 

The  distribution  of  q  and  »■  depends  on  k,  S^,  n^  and  r^. 
The  distribution  of  u*  depends  additionally  upon  p  and  S.  For 
fixed  choices  of  these  parameters  the  distributions  may  be 
determined  by  Monte  Carlo  sampling. 
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4.  INTERVAL  AND  MEDIAN  UNBIASED  POINT  ESTIMATES 

Two  sided  100  interval  estimates  of  6,  Y  and  xp(S) 

may  be  calculated  in  terms  of  the  ML  estimates  and  the  percentage 
points  of  q,  w*  and  u*  as: 


A  A 


‘^(l-a/2) 

<  8  <  ®\/2 

(10) 

*  '  "*l-a/2/8 

<  y  K  y  -  w*/2/S 

(11) 

xp(S)  •  exp  [-u^.^/g]  <  *p(S)  <  xp(S)  •  exp  [-u*a/2/S]  (12) 

Median  unbiased  point  estimates  are  calculable  as: 

A  A 


8'  *  8/V  50 

(13) 

~ 

V  •  Y  ■  w0. 50/8 

(14) 

x^(S)  -  xpCS)  •  exp  C-Uq_ S0/6J 

(IS) 

5.  PRECISION  MEASURES 

A  useful  index  of  how  precisely  0  is  determined  by  the  series 
of  life  tests  is  the  ratio  R  of  the  upper  to  lower  ends  of  its 
confidence  interval. 


R  “  q(l-a/2)/<l(a/2)  06) 

For  Y  the  median  length  LQ  of  the  confidence  interval  is 
a  convenient  index  calculated  as 

L0.S0  *  (wl-a/l  '  wa/2)/SqQ.SQ-  .  —  07) 
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For  Xp(S)  the  median  ratio  Rg  upper  to  lower  confidence 

interval  is  recommended. 

R0.S0  *  tC-“a/2  *  ul-a/2)/6VS01  U8) 

NUMERICAL  EXAMPLE 

Rolling  contact  fatigue  tests  of  n.^  =  10  hardened  steel 
specimens  were  conducted  at  k  3  4  levels  of  the  contact  stress. 

The  tests  were  continued  until  all  elements  had  failed  (r^  =  10). 
The  results  are: 


Stress 

106-psi 

Ordered  Lives 

0.87 

1.67,  2.20,  2.51,  3.0C,  3.90,  4.70,  7.53, 
27.76,  37.4 

14.70 

0.99 

0.80,  1.00,  1.37,  2.25,  2.95,  3.70,  6.07, 
7.05,  7.37 

6.65, 

1.09 

0.012,  0.18,  0.20,  0.24,  0.26,  0.32,  0.32 
0.44,  0.88 

,  0.42 

1.18 

0.073,  0.098,  0.117,  0.135,  0.175,  0.262, 
0.350,  0.386,  0.456 

0.270 

The  ML  estimates  are 

3  *  1.166  Y  *  13.89  n  *  2.20 

The  5-th,  50-th  and  95-th  percentiles  of  the  distribution 
of  q,  w*  and  u*  corresponding  to  k  ■  4,  n^  *  r.  *  10  and  the 
specified  values  of  S-  are  given  in  Table  1.  The  distribution 
of  u*  was  evaluated  for  p  *  0.10  with  S  3  0.75  x  10  psi  and 
with  the  four  test  stresses. 
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Percentiles  of  Pivotal  Functions 
k  -  4,  n  ■  r  -  10 
Sj-0.87,  S2-0.99,  Sj-1.09,  S4-1.18 


0.05 

0.50 

0.95 

q  -  B/8 

0.8459 

1.024 

1.277 

A  A 

w*  .  (y-y)0 

-2.433 

-0.3783 

2.293 

u*  ■  ®l0S  [x0.10/x0.10]; 

S  -  0.75 

-0.9238 

0.0555 

1.023 

S  -  S1  -  0.87 

-0.6495 

0.0441 

0.8209 

S  -  S2  -  0.99 

-0.S170 

0.0305 

0.7520 

S  -  S3  -  1.09 

-0.5309 

0.0318 

0.7671 

S  -  S4  -  1.18 

-0.6079 

0.0396 

0.8193 

90%  confidence  intervals  for  8  and  y  are 

0.913  -  1.166/1.277  <  S  <  1.166/0.8459  =  1.378 

11.92  -  13.889  -  2.293/1.166  <  y  <  13.889  +  2.433/1.166  *  1 

Median  unbiased  estimates  are 
S'  -  1.166/1.024  -  1.139 

y'  •  13.889  ♦  0.3783/1.166  -  14.21 


5.98 
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The  precision  measures  based  on  a  901  interval  are 
R  *  1.51 


SL0.5Q 


4.62 


The  confidence  intervals  and  unbiased  estimates  for  xQ  1Q 
are  listed  below  for  each  stress  along  with  the  precision  measure 


Stress 

Median  Unbiased 
x0  Estimate 

901  Confidence 
Lower  X°-10 

Interval 

Upper 

d8 

—0.50— 

0.75 

16.55 

7.22 

38.4 

6.70 

0.87 

2.13 

1.09 

3.86 

4.21 

0.99 

0.358 

0.193 

0.572 

3.42 

1.09 

0.094 

0.050 

0.152 

3.55 

1.18 

0.031 

0.016 

0.054 

4.03 

Figure  1 

shows  R^  plotted  against  stress 

.  The  minimum 

value  is  only 

slightly  larger 

then  for  a  single 

uncensored 

sample 

of  size  n  *  40. 
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